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Using data from the SUGAR and the AGASA experiments taken during a 10 yr period with nearly 
uniform exposure to the entire sky, we search for large-scale anisotropy patterns in the arrival 
directions of cosmic rays with energies > 10 19 ' 6 eV. We determine the angular power spectrum from 
an expansion in spherical harmonics for modes out to £ = 5. Based on available statistics, we find no 
significant deviation from isotropy. We compare the rather modest results which can be extracted 
from existing data samples with the results that should be forthcoming as new full-sky observatories 
begin operation. 



I. INTRODUCTION 

Ultra high energy cosmic rays are one of the most enig- 
matic phenomena in the universe. Despite the fact that 
the existence of particles with energies > 10 20 eV has 
been known for over 40 years, their origin continues to 
be an intriguing puzzle [1] . The distribution of arrival di- 
rections is perhaps the most helpful observable in yielding 
clues about cosmic ray origin. On the one hand, if cos- 
mic rays cluster within a small angular region [2] or show 
directional alignment with powerful compact objects [3], 
one might be able to associate them with isolated sources 
in the sky. On the other hand, if the distribution of 
arrival directions exhibits a large-scale anisotropy, this 
could indicate whether or not certain classes of sources 
are associated with large-scale structures (such as the 
Galactic plane or the Galactic halo). In this paper, we 
focus our attention on the search for such large-scale ce- 
lestial patterns. 

Cosmic ray air shower detectors which experience sta- 
ble operation over a period of a year or more will have 
a uniform exposure in right ascension, a. A traditional 
technique to search for large-scale anisotropies is then 
to fit the right ascension distribution of events to a sine 
wave with period 2n /m (m th harmonic) to determine the 
components (x, y) of the Rayleigh vector [4] 

2 N 2 N 

x = — cos(ma 4 ), y = — ^ sin(mai). (1) 

i=l i=l 

The TO th harmonic amplitude of TV measurements a, is 
given by the Rayleigh vector length 1Z = (x 2 + y 2 ) 1 ! 2 . 
The expected length of such vector for values ran- 
domly sampled from a uniform phase distribution is 
TZq = 2j\f~N. The chance probability of obtaining 
an amplitude with length larger than that measured is 
p{> Tty = e ~ k ° , where k = 1Z 2 /1Zq. To give a specific 
example, a vector of length fco > 6.6 would be required 
to claim an observation whose probability of arising from 
random fluctuation was 0.0013 (a "3er" result). For the 
ultra high energy (> 10 19 ' 6 eV) regime, all experiments 
to date have reported k <C 6.6 [5]. This does not imply 



an isotropic distribution, but it merely means that avail- 
able data are too sparse to claim a statistically significant 
measurement of anisotropy by any of these experiments. 
In other words, there may exist anisotropies at a level 
too low to discern given existing statistics [6] . 

The right harmonic analyses are completely blind to 
intensity variations which depend only on declination, 5. 
Combining anisotropy searches in a over a range of de- 
clinations could dilute the results, since significant but 
out of phase Rayleigh vectors from different declination 
bands can cancel each other out. Moreover, the analysis 
methods that consider distributions in one celestial coor- 
dinate, while integrating away the second, have proved to 
be potentially misleading [7]. An unambiguous interpre- 
tation of anisotropy data requires two ingredients: expo- 
sure to the full celestial sphere and analysis in terms of 
both celestial coordinates. Though the statistics are very 
limited at present, this article describes a first step in this 
direction. In the next section we combine data from the 
Sydney University Giant Air-shower Recorder (SUGAR) 
and the Akeno Giant Air Shower Array (AGASA) taken 
during a 10 yr period with nearly uniform exposure to 
the entire sky. After that, in Sec. Ill, we apply the power 
spectrum estimation technique [8] to interpret the distri- 
bution of arrival directions. Our conclusions are collected 
in Sec. IV. 



II. EXPERIMENTAL DATA SETS 

The SUGAR array was operated from January 1968 
to February 1979 in New South Wales (Australia) at a 
latitude of 30.5° South and longitude 149°38' East [9]. 
The array consisted of 47 independent stations on a rect- 
angular grid covering an area S ~ 70 km 2 . The primary 
energy was determined from the total number of muons, 
N^, traversing the detector at the measured zenith angle 
9. The total aperture for incident zenith angles between 
9i and 82 was found to be 

A= S p{N^6) cosOdfl. (2) 
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Here, p(Np, 9) is the probability that a shower falling 
within the physical area was detected, S cos 9 is the pro- 
jected surface of the array in the shower plane, and dQ is 
the acceptance solid angle. The SUGAR Collaboration 
reports [9] a reasonable accuracy in assessing the shower 
parameters up to 9 = 73°. The particulars of the events 
with primary energy > 10 19 6 eV are given in Table I. 
The estimated angular uncertainty for showers that trig- 
gered 5 or more stations is reported as 3° sec 9 [9] . How- 
ever, the majority of events were only viewed by 3 or 4 
stations, and for these the resolution appears to be as 
poor as 10° [10]. Of particular interest for this analysis, 
p(Nf, > 10 8 ,6» < 55°) « 0.85 [11], yielding a total aper- 
ture A w 125 km 2 sr. This provides an exposure rea- 
sonably matched to that of AGASA, which is described 
next. 

The AGASA experiment occupies farm land near the 
village of Akeno (Japan) at a longitude of 138°30' East 
and latitude 35°30' North [12]. The array, which con- 
sists of 111 surface detectors deployed over an area of 
about 100 km 2 , has been running since 1990. About 
95% of the surface detectors were operational from March 
to December 1991, and the array has been fully op- 
erational since then. A prototype detector operated 
from 1984 to 1990 and has been part of AGASA since 
1990 [13]. The aperture for events with primary zenith 
angle 0° < 9 < 45° and energies beyond 10 1925 eV is 
found to be A « 125 km 2 sr [12]. The angular resolu- 
tion for these events is 1.6° [14]. The arrival directions of 
cosmic rays with energy > 10 19 6 eV are given in Table II. 

The expected event rate is found to be 
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where (E 3 J(E)) w 10 24 6 cV 2 uT 2 s" 1 sr" 1 stands for 
the observed ultra high energy cosmic ray flux, which has 
a cutoff at E 2 = 10 20 5 cV [1]. Thus, in approximately 
10 yr of running each of these experiments should col- 
lect w 50 events above E\ = 10 19 6 eV, arriving with 
a zenith angle < 9 max . Here, 9 max = 45° for AGASA 
and # max = 55° for SUGAR. Our sub-sample for the full- 
sky anisotropy search consists of the 50 events detected 
by AGASA from May 1990 to May 2000 [6], and the 49 
events detected by SUGAR with 9 < 55° [11]. Note that 
we consider the full data sample for the 11 yr lifetime 
of SUGAR (in contrast to the 10 yr data sample from 
AGASA). This roughly compensates for the time vari- 
ation of the sensitive area of the experiment as detec- 
tors were deployed or inactivated for maintenance. The 
arrival directions of the 99 events are plotted in Fig. 1 
(equatorial coordinates B.1950). 

A detector at latitude ao that has continuous oper- 
ation with constant exposure in right ascension and is 
fully efficient for 9 < # max has relative exposure with the 
following dependence on declination [15] 
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FIG. 1. Arrival direction of the 99 events observed above 
10 19 ' 6 eV by the SUGAR (9 < 55°) and the AGASA (6 < 45°) 
experiments (equatorial coordinates B.1950). 



ui(5) cx (cosa cos<5 sin a max + a max sina sin<5) , (4) 

where a max , the local hour angle at which the zenith 
angle becomes equal to # max , is given by 



with 



if £ > 1 

7T if £ < -1 

cos -1 £ otherwise 
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We normalize the relative exposure by scaling the decli- 
nation angle dependence given in Eq. (4) by an estimate 
of the areas and operation times for SUGAR and AGASA 
(recall that the stated area of SUGAR is an overestimate, 
and this compensated to some degree by the longer run- 
ning time). The result, displayed in Fig. 2, shows that 
the combined exposure of these arrays is nearly uniform 
over the entire sky. 



III. CORRELATIONS AND POWER SPECTRUM 

We begin this section with a general introduction to 
the calculation of the angular power spectrum and the 
determination of the expected size of intensity fluctua- 
tions. The technique is then applied to the AGASA and 
SUGAR data in order to check for fluctuations beyond 
those expected from an isotropic distribution. 

Let us start by defining the directional phase space of 
the angular distribution of cosmic ray events in equato- 
rial coordinates, (a, S). (i) The direction of the event is 
described by a unit vector 

n = sin 5 (i cos a + j sin a) + k cos S ; (7) 

(ii) The solid angle is given by 
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FIG. 2. Declination dependence of SUGAR and AGASA 
relative exposures (dotted). The solid line indicates the com- 
bined declination dependence of the relative exposure. 



d 2 n = sin 5 d5 da ; (8) 

(iii) The delta function for the solid angle is defined as 

oo 

5(n,n') = 5{cos8 - cos 5') V] 5 (a - a' + , (9) 

m— — oo 

so that, as usual, 

/(n) = |5(n,n') /(n') d 2 n' ; (10) 

(iv) The probability distribution P(n)d 2 n of events can 
be employed for the purpose of computing the averages 



/ = / /(n) P(n) d 2 n ; 



(11) 



Finally, (v) for a sequence of N different cosmic ray 
events (ni, . . . , njv) one may assume an independent dis- 
tributions for each event, i.e. 



N 



N 



P N (m,..., n N ) I] d 2 n, = IJ{P(nO d 2 ^} . (12) 

i i 

For a sequence of events (ni, . . . , njv) let us describe 
the angular intensity as the random variable 



1 N 



(13) 



From Eqs. (12) and (13) it follows that 

7(H) = j ... j I(n)P N (n 1 ,... 1 n N )Y[d 2 n. 



= P(n) 



(14) 



The two point correlation function G(n, n') = 7(n)7(n / ) 
is defined via 



G( 



n,n') = j ...j J(n) 7(n') P w (m, . . . , n w ) J|d 2 n 



l<5(n,n')P(n)+[l-l) P(n)./>(n') 



(15) 

The "power spectrum" of the correlation function is 
determined by the eigenvalue equation 
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G(n,n') ^A(n') d 2 n' = A ^ A (n). 



(16) 



In this regard it is useful to introduce Dirac notation to 
indicate the inner product 

= I V» V(n) d 2 n • (17) 

With this in mind, Eq. (16) reads 

G |Va) - A |Va) • (18) 
In the limit of a large number of events N — > oo, 

lim G(n, n') = G M (n, n') = P(n)P(n') , (19) 

JV— >oo 



or equivalcntly, 



Goo = |P) (P| • 



(20) 



In such a limit, fluctuations can be neglected and we find 
only two possible values in the spectrum: (i) There is a 
non-degenerate non-zero eigenvalue 



with 



Goo |P) = Aoo |P) , 

Aoo - (P|P) = J P 2 (n)d 2 n. 



(21) 



(22) 



(ii) For every state |/) orthogonal to |P) with mean value 
/ = (P\f) = 0, there exists a zero eigenvalue in the power 
spectrum 

Goo \f) = |P) (P\f) = \ I Pfd 2 n\ |P) = f\P) = . 



(23) 

Let us now turn to consider the effects of finite N. Defin- 
ing the fluctuations in the intensity by 



A/(n) = /(n) - /(n) = 7(n) - P(n), (24) 
the two point correlation function can be re-written as 



G(n,n') = 7(n)7(n') = 7(n) J(n') + A/(n)A7(n') 



Goo(n,n') + A7(n)A7(n'), 



(25) 



with 
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FIG. 3. The angular power spectrum is indicated by the 
squares. The horizontal lines indicate the mean value, 
C(£) = (4nN)~ , expected for an isotropic distribution. The 
upper shaded band shows the la fluctuation around the mean 
value for N = 99. The band was obtained from 1000 sets of 
Monte Carlo simulations of 99 events each, including small 
corrections for uij . The resuts from the Monte Carlo are con- 
sistent with the errors computed using Eq. (48) for the case of 
uij = 1. For £ — 3, where there is a small excess compared to 
the expectation for isotropy, C(3) = 2.16 x 10~ 3 while the ex- 
pectation from a random distribution is Cmc(3) = 9.5 x 10 -4 , 
with a variance of 5.0 x 10~ 4 . The projected sensitivity for 
the Pierre Auger Observatory is also indicated on the plot by 
the lower shaded band. 
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G(n,n') Y lm (n')d 2 n' = \ tm Y lm (n) 



(29) 



where from Eq. (28) we have, 



_/(47r)- 1 if (*,m) = (0,0) ( . 

Alm ~ \ (AttN)- 1 if (i, m) £ (0, 0) ■ [M) 

The eigenf unctions form a useful set for expansions of the 
intensity over the celestial sphere 



J(n) = ^ ^2 ae m Y hn (n) . 



(31) 



e=o m=-l 



To incorporate the dependence on declination given in 
Eq. (4), let us re-define the angular intensity 



1 N 1 

7 ( n ) = >7 E-^ n ^). 



(32) 



where u>j is the relative exposure at arrival direction r\j 
given in Fig. 2 and Af is the sum of the weights ujJ 1 . 
Since the eigenvalues of the Yi m expansion are uniquely 
defined 



aim = J I(n) Yi m (n) d 2 n , 



(33) 



the replacement of Eq. (32) into Eq. (33) leads to the 
explicit form of the coefficients for our set of arrival di- 
rections 



at, i 
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(34) 



A/(n)A7(n') = - [6(n, n')P(n) - P(n)P(n') ] , (26) 

where Eq. (15) has been invoked. Putting all this to- 
gether, some general results follow: (i) For the N — > oo 
case, there is only one state with a finite eigenvalue 
Aoo, while the rest of the power spectrum corresponds 
to A = 0. (ii) For finite N, Eq. (26) implies that the 
fluctuations are of order A^ 1 . The power spectrum for 
large N then has one eigenvalue of order unity and the 
rest of the eigenvalues are of order A^ 1 . 
Now, for an isotropic distribution of n, 



(27) 



and the two point correlation function of Eq. (15) be- 
comes, 



1 

1 

N 



(28) 



The eigenvalue problem is solved by employing spherical 
harmonics [17] 



Squaring Eq. (34) and taking the average yields 

-g— O£ S m0 1 

a irn = —^ + 4^( 1 - S "> 5 "*> ■ 



(35) 



Equivalently, the mean square fluctuations of the coeffi- 
cients are determined by the power spectrum eigenvalues 
according to 



a 2 . = Xt 



(36) 



Although full anisotropy information is encoded into the 
coefficients ai m (tied to some specified coordinate sys- 
tem), the (coordinate independent) total power spectrum 
of fluctuations 



C{i) 



1 l 
+ 1) ^, 



(2£+l) 



(37) 



provides a gross summary of the features present in the 
celestial distribution together with the characteristic an- 
gular scale(s). Note that Eqs. (30) and (36) imply 
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if t = 
1 if^O 



(38) 



One must take care when combining data from exper- 
iments with differing angular resolutions, as in the case 
for SUGAR and AGASA. In this analysis, we consider 
only anisotropy on an angular scale larger than the the 
resolution of SUGAR, which has the worse resolution of 
the two. The power in mode £ is sensitive to variation 
over angular scales of £~ x radians [15]. Recalling that 
the estimated angular uncertainty for some of the events 
in the SUGAR sample is possibly as poor as 10° [10] we 
only look in this study for large scale patterns, going into 
the multipolc expansion out to £ — 5. 

Our results at this juncture are summarized in Fig. 3. 
The angular power spectrum is consistent with that ex- 
pected from a random distribution for all (analyzed) mul- 
tipoles, though there is a small (2a) excess in the data for 
I = 3. The majority of this excess comes from SUGAR 
data [18]. The decrease in error as £ increases may be un- 
derstood as a consequence of the fact that contributions 
to mode £ arise from variations over an angular scale £~ Y . 
If one compares to the expectation for isotropy, struc- 
tures characterized by a smaller angular scale, and hence 
larger £, can be ruled out with more significance than 
larger structures. 

To quantify the error, we study the fluctuations in C(£) 
for £ > 1. For simplicity, let us neglect the small effects 
of declination (viz., ojj = 1 Vj), and consider the random 
variable 



X, = 



C(£) 
C(£) 



A-kN 
2£+l 



E 



(39) 
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Denoting by P^(cos<5) the Legendre polynomial of order 
£ and employing the addition theorem for spherical har- 
monics, 



47T 



2£ + 



- J2 Y em (n)Y tm (n') = P e (n-n'), (40) 



Eqs. (34), (39), and (40) imply that 
2 



X, = 1 



N 



E p ^ ■ n j) • 



(41) 



l<i<j<N 



where Pg(l) = 1 has been invoked. Evidently, Xg = 1. 
Besides, 

= l + E E Pi(n t ■ n^nk ■ n q ) . 



l<i<j<N l<k<q<N 



(42) 



Since different pairs in the sum on the right hand side of 
Eq. (42) are uncorrelated for an isotropic distribution, it 
follows that 



(43) 



l<i<j<JV 



There are {N(N - l)/2} equivalent pairs in Eq. (43) 
which implies 



1 



X} = X t +2^1--)P £ ( ni .n 2 )2. (44) 
^From Eq. (40) we obtain 



P,(m • n 2 ) 2 = ( -i^-J £ y, m ( ni ) 2 Y lm (n 2 ) 
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Plugging Eq. (45) into Eq. (44) leads to 



1 \ 2 



AX " = {1 -NJ 2£+l 



or equivalcntly, 
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(45) 
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(47) 



yielding (for large N) 



/ ACf \ 2 
lim = — for £ > 1 , (48) 

N ^ \cf J 2£ + 1 ~ 

which is the variance on JQ. 



IV. OUTLOOK 

We have made a first full-sky anisotropy search using 
data from the SUGAR and AGASA experiments. At 
present, low statistics and poor angular resolution limits 
our ability to perform a very sensitive survey, but we can 
at least have a preliminary look at the first moments in 
the angular power spectrum. The data are consistent 
with isotropy, though there appears to be a small excess 
for £ = 3, arising mostly from the SUGAR data. 

There are two caveats in this analysis which should be 
kept in mind. First, from the published SUGAR results, 
it is difficult to make an exact determination of the expo- 
sure, as the sensitive area of the experiment varied as a 
function of time. Here, we assumed an area-time product 
of approximately 775 km 2 yr. Second, there is some un- 
certainty in the energy calibration. The SUGAR results 
are reported in terms of the number of vertical equivalent 
muons together with two possible models to convert this 
to primary energy. We have chosen the model yielding an 
energy spectrum which is in better agreement with the 
AGASA results [19]. It should be noted, though, that this 
spectrum does not agree well with the results of the High 
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Resolution Fly's Eye (HiRes) experiment [20]. Though 
there are uncertainties in the energy scale, the impact 
on this anisotropy search may not be so severe. This is 
because the energy cut of 10 19 6 eV is well above the last 
break in the spectral index at <~ 10 18 5 — 10 19 eV, and 
one would expect that all cosmic rays above this break 
share similar origins. 

In the near future we expect dramatically superior re- 
sults from the Pierre Auger Observatory. This observa- 
tory is designed to measure the energy and arrival direc- 
tion of ultra high energy cosmic rays with unprecedented 
precision. It will consist of two sites, one in the North- 
ern hemisphere and one in the Southern, each covering 
an area S ~ 3000 km 2 . The Southern site is currently 
under construction while the Northern site is pending. 
Once complete, these two sites together will provide the 
full sky coverage and well matched exposures which are 
crucial for anisotropy analyses. The base-line design of 
the detector includes a ground array consisting of 1600 
water Cerenkov detectors overlooked by 4 fluorescence 
eyes. The angular and energy resolutions of the ground 
arrays are typically less than 1.5° (multi-pole expansion 
t <~ 60) and less than 20%, respectively. The detec- 
tors are designed to be fully efficient (p w 1) out to 
#max = 60° beyond 10 19 eV, yielding a nearly uniform 
sky Am 1.4 x 10 4 km 2 sr [15]. In 10 yr of running the two 
arrays will collect « 4000 events above E\ — 10 19 6 eV. 
As can be seen from Fig. 3, such statistics will allow one 
to discern asymmetries at the level of about 1 in 10 4 . 
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TABLE I. Arrival directions of cosmic rays with energy > 10 19 ' 6 eV in equatorial coordinates (a, 8) together with the incident 
zenith angle (9) as reported by the SUGAR Collaboration [9]. 



9 [deg.] 


a [deg.] 


8 [deg.] 


9 [deg.] 


a [deg.] 


8 [deg.] 


9 [deg.] 


a [deg.] 


8 [deg.] 


70 


187.5 


31.8 


27 


333.0 


-56.6 


72 


354.0 


-75.2 


28 


117.0 


-3.3 


54 


231.0 


-31.0 


43 


357.0 


-57.4 


43 


147.0 


-43.6 


34 


130.5 


-27.4 


24 


121.5 


-32.0 


31 


105.0 


-38.7 


23 


288.0 


-51.4 


64 


204.0 


24.0 


29 


231.0 


-13.5 


28 


57.0 


-3.7 


35 


261.0 


-32.8 


38 


135.0 


4.3 


60 


340.5 


-43.2 


51 


300.0 


-41.9 


32 


238.5 


-9.5 


38 


16.5 


-68.5 


68 


181.5 


20.1 


31 


247.5 


-0.7 


33 


331.5 


-32.2 


23 


69 


-49.1 


40 


280.5 


-55.6 


64 


303 


13.6 


62 


126.0 


-39.9 


55 


226.5 


-31.1 


42 


316.5 


-65.1 


45 


189 


-15.5 


6 


114 


-25.8 


14 


145.5 


-38.1 


48 


154.5 


-27.7 


53 


268.5 


-81.7 


69 


216 


18.6 


46 


76.5 


10 


69 


315 


0.7 


67 


154.5 


3.1 


67 


123 


-38 


23 


277.5 


-8.5 


38 


160.5 


1.5 


53 


1.5 


-27.1 


54 


337.5 


-42.2 


57 


85.5 


20.5 


70 


88.5 


-64.6 


64 


303 


-47.2 


44 


142.5 


-24.6 


52 


172.5 


-30.8 


15 


99 


-24.7 


49 


30 


15.9 


49 


73.5 


0.2 


58 


355.5 


-0.5 


65 


283.5 


-50.6 


42 


19.5 


-62.8 


27 


331.5 


-15.5 


15 


144 


-33.1 


41 


12 


-14.1 


68 


163.5 


-5.9 


63 


27 


3.1 


46 


315 


-26.6 


56 


234 


1.2 


43 


94.5 


-41 


57 


229.5 


-78.7 


72 


117 


-35.9 


62 


262.5 


28.2 


22 


93 


-45.6 


59 


223.5 


-88.3 


40 


3 


-2.9 


33 


195 


-44.7 


59 


30 


-78.5 


38 


327 


-24.7 


23 


231 


-46.6 


52 


340.5 


-42.8 


41 


87 


6.1 


58 


183 


4.5 


57 


22.5 


22.2 


61 


244.5 


14.3 


65 


255 


-83 


57 


88.5 


-47.8 


57 


180 


-22.7 









TABLE II. AGASA events with 9 < 45° and mean energy > 10 196 eV in equatorial coordinates (a, 8) based on the B.1950 
reference frame. The first 2 rows refer to events measured using the prototype detector before 1990 [16]. 



a [deg.] 


8 [deg.] 


a [deg.] 


8 [deg.] 


a [deg.] 


8 [deg.] 


a [deg.] 


8 [deg. 


334.70 


38.15 


276.80 


35.26 


68.71 


30.00 


210.03 


50.14 


328.68 


27.36 


206.45 


34.95 


86.65 


58.48 






243.58 


-7.08 


236.31 


41.15 


284.40 


47.73 


53.01 


69.33 


287.03 


77.12 


2.29 


78.32 


143.23 


38.82 


267.68 


47.91 


7.85 


17.42 


255.27 


31.47 


171.55 


57.37 


100.17 


34.95 


123.54 


16.95 


208.34 


60.04 


28.28 


28.75 


18.25 


49.74 


18.07 


20.83 


280.91 


48.25 


73.27 


17.92 


167.77 


57.87 


263.60 


-1.57 


192.39 


30.87 


17.83 


19.73 


69.46 


29.80 


240.96 


23.13 


57.24 


26.95 


269.34 


74.10 


198.98 


53.16 


297.94 


18.57 


323.63 


7.87 


247.29 


34.70 


213.73 


37.93 


294.35 


70.98 


293.83 


-5.91 


33.82 


13.57 


166.55 


42.07 


146.04 


23.93 


258.77 


56.35 


167.52 


56.27 


348.37 


12.03 


293.67 


50.59 


55.38 


44.74 


287.13 


5.22 


113.95 


32.31 


55.59 


49.34 


345.16 


33.63 


339.45 


42.34 


68.59 


5.00 


59.55 


51.56 


211.47 


37.34 
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